#include  <stdio.h>
#include  <math.h>
#include  <stdlib.h>
#include <sys/io.h>
#include <unistd.h>
#include <string.h>
#include <stdbool.h>
//#include <syslib.h>
#define Npix  6
#define Notpix 50  //10 arcsec
#define Npixref  5
main()
{
    int     i,j,k,m,N,N1,N2,N3,N4,p;
    float	rxc1,ryc2,rRA3,rDEC4,rB15,rB26,rR17,rR28,rI9;
    float	cxc1,cyc2,cRA3,cDEC4,cB15,cB26,cR17,cR28,cI9;
    float   xc1,yc2,mag3,merr4,elli5,class6,fwhm7;
    float   otxc1,otyc2,otRA3,otDEC4,magB2,magR2,magI;
    int     Obj_label;
    char    ottype5[200];
    int     flag8;
	float   delx,dely,delr,delotx,deloty,delotr;
//    float   sum_delmag_R2,sum_delmag_B2,average_Deltamag_R2,average_Deltamag_B2;
    float   sum_B2,sum_R2,sum_comp,sum_I;
    float   average_comp,average_R2,average_B2,average_I;
    float   x[1000], y[1000],mag[1000],merr[1000],elli[1000],class[1000],fwhm[1000];
    int     flag[1000],otObjLabel[1000];
    float   ref1[1000],ref2[1000],ref3[1000],ref4[1000],ref5[1000],ref6[1000],ref7[1000],ref8[1000],ref9[1000];
    float   comp1[1000],comp2[1000],comp3[1000],comp4[1000],comp5[1000],comp6[1000],comp7[1000],comp8[1000],comp9[1000];
    float   otx[100],oty[100],otRA[100],otDEC[100],otmagB[100],otmagR[100],otmagI[100];
    char    ottype[100][200];
    FILE    *fp1,*fnomatch,*fcomp,*fdeltamag;
	char    *refall,*refcomp, *newimg,*otcoord;

	refall="refusno.cat";
	refcomp="refcomp.cat";
	newimg="newimage.cat";
    otcoord="otcoord.cat";

	//printf("refall = %s\n",refall);
	//printf("newimg = %s\n",newimg);
    //printf("otcoord =%s\n",otcoord);

	i=j=m=N=k=N1=N2=N3=N4=p=0;
    sum_comp=sum_R2=sum_B2=sum_I=0.0;
    fcomp=fopen("comp_c.db","w+");
	fnomatch=fopen("nomatch_c.db","w+");
    fdeltamag=fopen("deltamag_c.db","w+");
	fp1=fopen(newimg,"r");	
	if(fp1)
	{
	        while((fscanf(fp1,"%f %f %f %f %f %f %f %d\n",&xc1,&yc2,&mag3,&merr4,&elli5, &class6, &fwhm7, &flag8))!=EOF)
	        {
                x[i]=xc1;
                y[i]=yc2;
                mag[i]=mag3;
                merr[i]=merr4;
                elli[i]=elli5;
                class[i]=class6;
                fwhm[i]=fwhm7;
                flag[i]=flag8;
		        i++;
		    }
		N1=i;
	}
	fclose(fp1);
	
	fp1=fopen(otcoord,"r");	
	if(fp1)
	{
	        while((fscanf(fp1,"%f %f %f %f %f %f %f %s %d\n",&otxc1,&otyc2,&otRA3,&otDEC4,&magB2,&magR2,&magI, ottype5, &Obj_label))!=EOF)
	        {
                otx[m]=otxc1;
                oty[m]=otyc2;
                otRA[m]=otRA3;
                otDEC[m]=otDEC4;
                otmagB[m]=magB2;
                otmagR[m]=magR2;
                otmagI[m]=magI;
		        sprintf(ottype[m],"%s",ottype5);
                otObjLabel[m]=Obj_label;
                m++;
            }
        N3=m;
	}
	fclose(fp1);
	
//	 printf("####\n");
    //	 refall is for the usno catalog
	fp1=fopen(refall,"r");
        if(fp1)
        {
          while((fscanf(fp1,"%f %f %f %f %f %f %f %f %f\n",&rxc1,&ryc2,&rRA3,&rDEC4,&rB15,&rB26,&rR17,&rR28,&rI9))!=EOF)
          {
          ref1[k]=rxc1;
          ref2[k]=ryc2;
          ref3[k]=rRA3;
	      ref4[k]=rDEC4;	
	      ref5[k]=rB15;
	      ref6[k]=rB26;
	      ref7[k]=rR17;
	      ref8[k]=rR28;
          ref9[k]=rI9;
//          printf("%.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f\n",ref1[j],ref2[j],ref3[j],ref4[j],ref5[j],ref6[j],ref7[j],ref8[j],ref9[j]);
          k++;
          }
          N4=k;
        }
        fclose(fp1);

        
       	fp1=fopen(refcomp,"r");
        if(fp1)
        {
          while((fscanf(fp1,"%f %f %f %f %f %f %f %f %f\n",&cxc1,&cyc2,&cRA3,&cDEC4,&cB15,&cB26,&cR17,&cR28,&cI9))!=EOF)
          {
          comp1[j]=cxc1;
          comp2[j]=cyc2;
          comp3[j]=cRA3;
	      comp4[j]=cDEC4;	
	      comp5[j]=cB15;
	      comp6[j]=cB26;
	      comp7[j]=cR17;
	      comp8[j]=cR28;
          comp9[j]=cI9;
          j++;
          }
          N2=j;
        }
        fclose(fp1);




    for(m=0;m<N3;m++) //OT 
    {
	    for(i=0;i<N1;i++)  //newimg
	    {
            delotx=x[i]-otx[m];
            deloty=y[i]-oty[m];
	        delotr=sqrt(pow(delotx,2)+pow(deloty,2));
            p=0;
            if(delotr<Notpix)
            {
            //    for(k=0;k<N4;k++)
            //    {
            //        delx=x[i]-ref1[k];
            //        dely=y[i]-ref2[k];
	        //        delr=sqrt(pow(delx,2)+pow(dely,2));
            //        if(delr>Npixref) //this obj is in the usno and in the error circular of miniOT
            //        {
            //            p++;
            //     //       printf("x=%.3f y=%.3f ref1=%.3f ref2=%.3f !!!!!!!!!!!111delr=%.3f\n",x[i],y[i],ref1[k],ref2[k],delr);
            //        }
            //        
            //    }
            //   // printf("=== p=%d  N4=%d\n",p,N4);
            //    if(p==N4)
            //    {
	    	        fprintf(fnomatch,"%.3f %.3f %.3f %.3f %.3f %.3f %.3f %d %.3f %.3f %.3f %s %d\n",x[i],y[i],mag[i],merr[i],elli[i],class[i],fwhm[i],flag[i],otmagB[m],otmagR[m],otmagI[m],ottype[m],otObjLabel[m]);	
	        //    }	
            }    

        }
    }    
	
    for(j=0;j<N2;j++)
    {
	    for(i=0;i<N1;i++)
		{
			delx=x[i]-comp1[j];
			dely=y[i]-comp2[j];
			delr=sqrt(pow(delx,2)+pow(dely,2));
            if(delr<Npix)  //the obj in new image matched an obj in the compar obj list, comparason stars
			{
                    N++;
                    sum_comp=sum_comp+mag[i];
                    sum_R2=sum_R2+comp8[j];
                    sum_B2=sum_B2+comp6[j];
                    sum_I=sum_I+comp9[j];

				    fprintf(fcomp,"%.3f %.3f %.3f %.3f %.3f %.3f %.3f %d %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f\n",x[i],y[i],mag[i],merr[i],elli[i],class[i],fwhm[i],flag[i],comp1[j],comp2[j],comp3[j],comp4[j],comp5[j],comp6[j],comp7[j],comp8[j],comp9[j],delr);
			}
		}
    }
    
    if(N>0)  //have comp stars
    {
        average_comp=sum_comp/N;
        average_R2=sum_R2/N;
        average_B2=sum_B2/N;
        average_I=sum_I/N;
	    fprintf(fdeltamag,"%.3f  %.3f %.3f %.3f %d\n",average_comp, average_B2, average_R2, average_I, N);	
    }
    else  //no comp stars
    {
        average_comp=average_B2=average_R2=average_I=0.0;
	    fprintf(fdeltamag,"%.3f  %.3f %.3f %.3f %d\n",average_comp, average_B2, average_R2, average_I, N);	
    }
    
	fclose(fnomatch);
    fclose(fcomp);
    fclose(fdeltamag);
}
